#ifndef _INTERPOLATION_H_
#define _INTERPOLATION_H_

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <cstdlib>

#define maximal 1e10//无穷大
const static double epsl = 1e-8;//精度

int factorial(const int &k){//阶乘
    int p = 1;
    for (int i = 1; i <= p; i++){
        p = p*i;
    }
    return p;
}

class Function{//虚拟函数类
public:
    virtual double operator () (const double &x) const = 0;
    virtual double dev( const double &x) const{
        return ((*this)(x+epsl)-(*this)(x-epsl))/(2*epsl);}
    virtual double dev( const double &x, const int &n) const{//n阶导
        if( n == 0) return (*this)(x);
        else if( n > 0) return ((*this).dev(x+epsl,n-1)-(*this).dev(x-epsl,n-1))/(2*epsl);
        else {
            std::cerr << "Invalid input, meaningless return" << std::endl;
            return maximal;}
    } 
};

class zero : public Function{//零函数，用于测试
public:
    double operator () (const double &x) const {return 0;}
};



class Polynomial{
protected:
    std::vector<double> x;
    std::vector<double> fvalue;
    std::vector<std::vector<double>> coef;
    int n = 0;
public:
    Polynomial( Function &f, std::vector<double> x):
    x(x), n(x.size()) {
        int cnt = 0;
        for( int i = 0; i < n; i++){
            fvalue.push_back(f.dev(x[i],cnt));
            if( x[i] == x[i+1]) cnt++;
            else cnt = 0;
        }
    };
    Polynomial( std::vector<double> x, std::vector<double> fx):
    x(x), n(x.size()), fvalue(fx) {
        if(!check()){
            std::cerr << "Invalid input!" << std::endl;
            exit(0);
        }
    };
    bool check(){
        return (n == fvalue.size() && n > 0);
    }
    virtual void poly() = 0;
};

class Hermite : public Polynomial{
public:
    Hermite( Function &f, std::vector<double> x): Polynomial(f,x) {};
    Hermite( std::vector<double> x, std::vector<double> fx):
    Polynomial(x,fx) {};
    void poly(){
        std::vector<double> tmp;
        tmp.push_back(fvalue[0]);
        coef.push_back(tmp);
        for ( int i = 1, cnt = 0; i < n; i++)
        {
            if (x[i] == x[i - 1]) cnt++;
            else cnt = 0;
            tmp.clear();
            for (int k = 0; k <= cnt; k++)
            {
                tmp.push_back(fvalue[i - cnt + k] / factorial(k));
            }
            for (int j = cnt + 1; j <= i; j++)
            {
                double deva = (tmp[j-1]-coef[i-1][j-1])/(x[i]-x[i-j]);
                tmp.push_back(deva);
            }
            coef.push_back(tmp);
        }
    }
    std::vector<double> vectorprint(){
        std::vector<double> tmp;
        for( int i = 0; i < n; i++){
            tmp.push_back(coef[i][i]);
        }
        return tmp;
    }
    void print(){
        for (int i = 0; i < n; i++){
            std::cout << std::setw(8) << std::setiosflags(std::ios::left) << std::setprecision(2) << coef[i][i];
        }
        std::cout << std::endl;
    }
    void printtable(){
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j <= i; j++)
            {
                std::cout << std::setw(8) << std::setiosflags(std::ios::left) << coef[i][j] ;
            }
            std::cout << std::endl;
        }
    }
    void printpoly(){
        std::cout.precision(2);
        std::cout << coef[0][0];//要改起来实在麻烦，这里第一位就算是0也强制保留
        for( int i = 1; i < n; i++){
            if( coef[i][i] > 0){
                if( coef[i][i] > 0.001)
                    std::cout << "+" << coef[i][i];
                else{
                    std::cout << "+" << std::setiosflags(std::ios::scientific) << coef[i][i];
                    std::cout.unsetf(std::ios::scientific);
                    }
                }
            else if( coef[i][i] < 0){
                if( coef[i][i] < -0.001)
                std::cout << coef[i][i];
                else{
                    std::cout << std::setiosflags(std::ios::scientific) << coef[i][i];
                    std::cout.unsetf(std::ios::scientific);
                }
            }
            if( coef[i][i])    
                for( int j = 0; j < i; j++)
                {
                    if( x[j] > 0)
                        std::cout << "(x"  << -x[j] << ")";
                    else if( x[j] < 0)
                        std::cout << "(x" << "+"  << -x[j] << ")";
                    else
                        std::cout << "*x";
                }
        }
        std::cout << std::endl;
    }
    void printx(){
        for(auto a : x){
            std::cout << std::setw(8) << std::setiosflags(std::ios::left) << a;
        }
        std::cout << std::endl;
    }
    void printf(){
        for(auto a : fvalue){
            std::cout << std::setw(8) << std::setiosflags(std::ios::left) << std::setprecision(2) << a;
        }
        std::cout << std::endl;
    }
};


#endif